Background

The single-cell dataset generated by Bassez et al. (2021) Nature Medicine contains single-cell data from breast carcinoma biopsies from 29 patients taken immediately before (pre-treatment) and 9 +/- 2 days after (on-treatment) receiving anti-PD1 immunotherapy. Lacking an objective measure of clinical benefit following therapy, in this study the authors use T cell expansion following treatment as a proxy for response.

Full Data Original Source

Case Study Dataset Internal link

Goal

To compare the T cell subtype composition between responders and non-responders using ProjecTILs and reference maps for CD8 T cells and CD4 T cells.

R Environment

library(renv)
renv::restore()

library(patchwork)
library(ggplot2)
library(reshape2)
library(Seurat)
library(ProjecTILs)

scRNA-seq data preparation

ddir <- "input/Bassez_breast_data"
if (!file.exists(ddir)) {
    dir.create(ddir)
    dataUrl <- "https://drive.switch.ch/index.php/s/7jdqnvPZuJrriZB/download"
    download.file(dataUrl, paste0(ddir, "/tmp.zip"))
    unzip(paste0(ddir, "/tmp.zip"), exdir = ddir)
    file.remove(paste0(ddir, "/tmp.zip"))
}
# Count matrices
f1 <- sprintf("%s/1864-counts_tcell_cohort1.rds", ddir)

cohort1 <- readRDS(f1)
dim(cohort1)

[1] 25288 53382

# Metadata
meta1 <- read.csv(sprintf("%s/1870-BIOKEY_metaData_tcells_cohort1_web.csv", ddir))
rownames(meta1) <- meta1$Cell

head(meta1)

                                                             Cell nCount_RNA
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1 BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1       1252
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1 BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1       1869
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1 BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1       1000
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1 BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1       1288
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1 BIOKEY_13_Pre_AAAGATGGTTTGACAC-1       2056
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1 BIOKEY_13_Pre_AAAGATGTCAACGCTA-1       1224
                                 nFeature_RNA patient_id timepoint expansion
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1          700  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1          995  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1          627  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1          681  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1          789  BIOKEY_13       Pre       n/a
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1          634  BIOKEY_13       Pre       n/a
                                 BC_type cellSubType          cohort
BIOKEY_13_Pre_AAACCTGCAAGAAGAG-1   HER2+         gdT treatment_naive
BIOKEY_13_Pre_AAACGGGTCGTAGGAG-1   HER2+     NK_REST treatment_naive
BIOKEY_13_Pre_AAAGATGAGCAGGCTA-1   HER2+       CD4_N treatment_naive
BIOKEY_13_Pre_AAAGATGCAAGCCGTC-1   HER2+      CD8_EM treatment_naive
BIOKEY_13_Pre_AAAGATGGTTTGACAC-1   HER2+       CD8_N treatment_naive
BIOKEY_13_Pre_AAAGATGTCAACGCTA-1   HER2+      CD8_RM treatment_naive

data.seurat <- CreateSeuratObject(cohort1, project = "Cohort1_IT", meta.data = meta1)

Subset pre-treatment biopsies

table(data.seurat$timepoint)

   On   Pre 
27528 25854 

data.seurat <- subset(data.seurat, subset = timepoint == "Pre")

Subset T cells according to annotation by the authors

table(data.seurat$cellSubType)

               CD4_EM                CD4_EX  CD4_EX_Proliferating 
                 4810                  1707                   107 
                CD4_N               CD4_REG CD4_REG_Proliferating 
                 3007                  2861                    69 
               CD8_EM              CD8_EMRA                CD8_EX 
                 5443                   290                  2286 
 CD8_EX_Proliferating                 CD8_N                CD8_RM 
                  340                   354                  1963 
                  gdT               NK_CYTO               NK_REST 
                 1000                   230                   938 
           Vg9Vd2_gdT 
                  449 

data.seurat <- subset(data.seurat, subset = cellSubType %in% c("NK_REST", "Vg9Vd2_gdT",
    "gdT", "NK_CYTO"), invert = T)

We will remove samples that are too small and downsample large ones, in order to have similar contribution from all patients/samples. Downsampling large samples also speeds up downstream computations.

ds <- 1000
min.cells <- 200

tab <- table(data.seurat$patient_id)
keep <- names(tab)[tab > min.cells]
data.seurat <- subset(data.seurat, patient_id %in% keep)


Idents(data.seurat) <- "patient_id"
data.seurat <- subset(data.seurat, cells = WhichCells(data.seurat, downsample = ds))
table(data.seurat$patient_id)

 BIOKEY_1 BIOKEY_10 BIOKEY_11 BIOKEY_12 BIOKEY_13 BIOKEY_14 BIOKEY_15 BIOKEY_16 
     1000      1000       555      1000      1000       939      1000      1000 
BIOKEY_19  BIOKEY_2 BIOKEY_21 BIOKEY_24 BIOKEY_27 BIOKEY_28  BIOKEY_3 BIOKEY_31 
     1000       960       203       347       601       594       266       349 
 BIOKEY_4  BIOKEY_5  BIOKEY_6 
     1000      1000       616 

ProjecTILs analysis

options(timeout = max(900, getOption("timeout")))

download.file("https://figshare.com/ndownloader/files/38921366", destfile = "CD8T_human_ref_v1.rds")
ref.cd8 <- load.reference.map("CD8T_human_ref_v1.rds")

[1] "Loading Custom Reference Atlas..."
[1] "Loaded Custom Reference map Human CD8 TILs"

download.file("https://figshare.com/ndownloader/files/39012395", destfile = "CD4T_human_ref_v1.rds")
ref.cd4 <- load.reference.map("CD4T_human_ref_v1.rds")

[1] "Loading Custom Reference Atlas..."
[1] "Loaded Custom Reference map custom_reference"

a <- DimPlot(ref.cd8, cols = ref.cd8@misc$atlas.palette, label = T) + theme(aspect.ratio = 1) +
    ggtitle("CD8 T reference") + NoLegend()

b <- DimPlot(ref.cd4, cols = ref.cd4@misc$atlas.palette, label = T) + theme(aspect.ratio = 1) +
    ggtitle("CD4 T reference") + NoLegend()

a | b

ncores = 8
cd4.projected <- Run.ProjecTILs(data.seurat, ref.cd4, ncores = ncores, split.by = "patient_id")

cd8.projected <- Run.ProjecTILs(data.seurat, ref.cd8, ncores = ncores, split.by = "patient_id")

s1 <- subset(cd4.projected, subset = patient_id == "BIOKEY_13")
a <- plot.projection(ref.cd4, query = s1, linesize = 0.5, pointsize = 0.2) + ggtitle("Patient 13 - CD4 T")

s2 <- subset(cd8.projected, subset = patient_id == "BIOKEY_13")
b <- plot.projection(ref.cd8, query = s2, linesize = 0.5, pointsize = 0.2) + ggtitle("Patient 13 - CD8 T")

s3 <- subset(cd4.projected, subset = patient_id == "BIOKEY_10")
c <- plot.projection(ref.cd4, query = s3, linesize = 0.5, pointsize = 0.2) + ggtitle("Patient 10 - CD4 T")

s4 <- subset(cd8.projected, subset = patient_id == "BIOKEY_10")
d <- plot.projection(ref.cd8, query = s4, linesize = 0.5, pointsize = 0.2) + ggtitle("Patient 10 - CD8 T")

(a | b)/(c | d)

genes4radar <- c("CD4", "CD8A", "TCF7", "CCR7", "IL7R", "LMNA", "GZMA", "GZMK", "FGFBP2",
    "XCL1", "CRTAM", "TOX", "PDCD1", "HAVCR2", "PRF1", "GLNY", "GZMB", "TRAV1-2",
    "KLRB1", "FOXP3")


plot.states.radar(ref = ref.cd8, query = cd8.projected, genes4radar = genes4radar)

plot.states.radar(ref = ref.cd4, query = cd4.projected, genes4radar = genes4radar)

data.seurat$functional.cluster.cd4 <- NA
data.seurat@meta.data[colnames(cd4.projected), "functional.cluster.cd4"] <- cd4.projected$functional.cluster

data.seurat$functional.cluster.cd8 <- NA
data.seurat@meta.data[colnames(cd8.projected), "functional.cluster.cd8"] <- cd8.projected$functional.cluster

annos <- data.seurat@meta.data[, c("functional.cluster.cd4", "functional.cluster.cd8")]

consensus <- apply(annos, 1, function(x) {
    if (is.na(x[["functional.cluster.cd4"]]) & !is.na(x[["functional.cluster.cd8"]])) {
        x[["functional.cluster.cd8"]]
    } else if (is.na(x[["functional.cluster.cd8"]]) & !is.na(x[["functional.cluster.cd4"]])) {
        x[["functional.cluster.cd4"]]
    } else {
        NA
    }
})

data.seurat$functional.cluster <- consensus
data.seurat$functional.cluster <- factor(data.seurat$functional.cluster, levels = unique(data.seurat$functional.cluster))

a <- ref.cd8@misc$atlas.palette
b <- ref.cd4@misc$atlas.palette
cols <- c(a, b)

by.exp <- SplitObject(data.seurat, split.by = "expansion")

tb <- table(factor(by.exp$E$functional.cluster, levels = names(cols)))
tb <- tb * 100/sum(tb)
tb.m <- reshape2::melt(tb)
colnames(tb.m) <- c("Cell_state", "Perc_cells")

p1 <- ggplot(tb.m, aes(x = Cell_state, y = Perc_cells, fill = Cell_state)) + geom_bar(stat = "identity") +
    theme_bw() + scale_fill_manual(values = cols) + theme(axis.text.x = element_blank(),
    legend.position = "none") + ggtitle("Responders") + ylim(0, 40)


tb <- table(factor(by.exp$NE$functional.cluster, levels = names(cols)))
tb <- tb * 100/sum(tb)
tb.m <- reshape2::melt(tb)
colnames(tb.m) <- c("Cell_state", "Perc_cells")

p2 <- ggplot(tb.m, aes(x = Cell_state, y = Perc_cells, fill = Cell_state)) + geom_bar(stat = "identity") +
    theme_bw() + scale_fill_manual(values = cols) + theme(axis.text.x = element_blank(),
    legend.position = "right") + ggtitle("Non-Responders") + ylim(0, 40)

p1 | p2

which.types <- table(data.seurat$functional.cluster) > 100  # disregard small populations

states_all <- names(cols)

cols_use <- cols[names(which.types)][which.types]

norm.c <- table(by.exp$NE$functional.cluster)/sum(table(by.exp$NE$functional.cluster))
norm.q <- table(by.exp$E$functional.cluster)/sum(table(by.exp$E$functional.cluster))

foldchange <- norm.q[which.types]/norm.c[which.types]
foldchange <- sort(foldchange, decreasing = T)

tb.m <- reshape2::melt(foldchange)
colnames(tb.m) <- c("Cell_state", "Fold_change")
pll <- list()
ggplot(tb.m, aes(x = Cell_state, y = Fold_change, fill = Cell_state)) + geom_bar(stat = "identity") +
    scale_fill_manual(values = cols_use) + theme_bw() + geom_hline(yintercept = 1) +
    scale_y_continuous(trans = "log2") + theme(axis.text.x = element_blank(), legend.position = "left") +
    ggtitle("Responders vs. Non-responders")

Correspondence analysis of subtype composition

Do patients cluster based on their subtype composition? are subtype composition “classes” related to response to ICB? we can perform correspondence analysis (CA) to try to answer these questions, and group patients based on their composition of T cell subtypes.

library("FactoMineR")
min.cells <- 30

# Only use samples with sufficient amount of cells
tab <- table(data.seurat$patient_id)
samples.use <- names(tab)[tab > min.cells]
data.seurat <- subset(data.seurat, subset = patient_id %in% samples.use)

subtype_comp_table <- prop.table(table(data.seurat$functional.cluster, data.seurat$patient_id),
    margin = 2)
res.CA <- CA(t(subtype_comp_table), graph = FALSE)
a <- plot(res.CA) + theme(aspect.ratio = 1)

# Color by response
coords <- as.data.frame(res.CA$row$coord)
coords.ct <- as.data.frame(res.CA$col$coord)
coords.all <- rbind(coords[, 1:5], coords.ct[, 1:5])

is.R <- unique(data.seurat$patient_id[data.seurat$expansion == "E"])
is.NR <- unique(data.seurat$patient_id[data.seurat$expansion == "NE"])
coords[["Responder"]] <- "n/a"
coords$Responder[rownames(coords) %in% is.R] <- "Responder"
coords$Responder[rownames(coords) %in% is.NR] <- "Non-responder"
coords[["Patient"]] <- rownames(coords)

b <- ggplot(coords, aes(x = `Dim 1`, y = `Dim 2`, color = Responder)) + geom_point(size = 4) +
    theme_bw() + theme(aspect.ratio = 1) + scale_color_manual(values = c("#E9E9E9",
    "#fb0505", "#08e579")) + xlim(min(coords.all[, "Dim 1"]), max(coords.all[, "Dim 1"])) +
    ylim(min(coords.all[, "Dim 2"]), max(coords.all[, "Dim 2"]))

a | b

Discussion

Our automated ProjecTILs analysis of this patient cohort indicate that compared to non-responders, early responders to PD-1 blockade (in terms of T cell expansion) have higher relative amounts of exhausted T cells (both CD8+ and CD4+) infiltrating their tumors, in agreement with previous studies (Daud et al. 2016; Thommen et al. 2018; Kumagai et al. 2020; Bassez et al. 2020). This supports the notion that patients with a pre-existing tumor-specific CD8 T cell response in their tumors are more likely to respond to PD-1 blockade therapy.

Because tumoricidal CD8 Tex are derived from Tpex (Siddiqui et al. 2019; Miller et al. 2019), the presence of CD8 Tex seems to be a good proxy for the presence of smaller populations of quiescent Tpex, that have the ability to proliferate and differentiate into Tex cells following PD-1 blockade (Siddiqui et al. 2019; Miller et al. 2019). Indeed we observed a small population of Tpex-like cells, with increased levels R vs RN. Because there are very few Tpex their gene profile might be less robust than others, but compared to Tex they display higher levels of Tpex-associated genes, such as XCL1 and CD200, and lower levels of terminal-exhaustion markers such as HAVCR2 and ENTPD1.

In their study, Bassez et al. also observed increased levels in R vs NR of a population they referred to as exhausted CD4 T cells. Here we observed also a trend for increased levels of follicular-helper CD4 T cells, that display exhaustion features (e.g. TOX, PDCD1, ENTPD1 expression), although the most prominent enrichment in ICB-responders is for exhausted CD8 T cells.

Further reading

Dataset original publication - Bassez et al

ProjecTILs case studies - INDEX - Repository

The ProjecTILs method Andreatta et. al (2021) Nat. Comm. and code

References

• Bassez, A., Vos, H., Van Dyck, L. et al. A single-cell map of intratumoral changes during anti-PD1 treatment of patients with breast cancer. Nat Med 27, 820–832 (2021). https://doi.org/10.1038/s41591-021-01323-8

• Daud, A.I., Loo, K., Pauli, M.L., Sanchez-Rodriguez, R., Sandoval, P.M., Taravati, K., Tsai, K., Nosrati, A., Nardo, L., Alvarado, M.D., Algazi, A.P., Pampaloni, M.H., Lobach, I.V., Hwang, J., Pierce, R.H., Gratz, I.K., Krummel, M.F., Rosenblum, M.D., 2016. Tumor immune profiling predicts response to anti-PD-1 therapy in human melanoma. J. Clin. Invest. 126, 3447–3452. https://doi.org/10.1172/JCI87324

• Gros, A., Robbins, P.F., Yao, X., Li, Y.F., Turcotte, S., Tran, E., Wunderlich, J.R., Mixon, A., Farid, S., Dudley, M.E., Hanada, K., Almeida, J.R., Darko, S., Douek, D.C., Yang, J.C., Rosenberg, S. a, 2014. PD-1 identifies the patient-specific in filtrating human tumors. J. Clin. Invest. 124, 2246–59. https://doi.org/10.1172/JCI73639.2246

• Miller, B.C., Sen, D.R., Al Abosy, R., Bi, K., Virkud, Y.V., LaFleur, M.W., Yates, K.B., Lako, A., Felt, K., Naik, G.S., Manos, M., Gjini, E., Kuchroo, J.R., Ishizuka, J.J., Collier, J.L., Griffin, G.K., Maleri, S., Comstock, D.E., Weiss, S.A., Brown, F.D., Panda, A., Zimmer, M.D., Manguso, R.T., Hodi, F.S., Rodig, S.J., Sharpe, A.H., Haining, W.N., 2019. Subsets of exhausted CD8+ T cells differentially mediate tumor control and respond to checkpoint blockade. Nat. Immunol. 20, 326–336. https://doi.org/10.1038/s41590-019-0312-6

• Siddiqui, I., Schaeuble, K., Chennupati, V., Fuertes Marraco, S.A., Calderon-Copete, S., Pais Ferreira, D., Carmona, S.J., Scarpellino, L., Gfeller, D., Pradervand, S., Luther, S.A., Speiser, D.E., Held, W., 2019. Intratumoral Tcf1+PD-1+CD8+ T Cells with Stem-like Properties Promote Tumor Control in Response to Vaccination and Checkpoint Blockade Immunotherapy. Immunity 50, 195–211.e10. https://doi.org/10.1016/J.IMMUNI.2018.12.021

• Thommen, D.S., Koelzer, V.H., Herzig, P., Roller, A., Trefny, M., Dimeloe, S., Kiialainen, A., Hanhart, J., Schill, C., Hess, C., Prince, S.S., Wiese, M., Lardinois, D., Ho, P.C., Klein, C., Karanikas, V., Mertz, K.D., Schumacher, T.N., Zippelius, A., 2018. A transcriptionally and functionally distinct pd-1 + cd8 + t cell pool with predictive potential in non-small-cell lung cancer treated with pd-1 blockade. Nat. Med. 24, 994. https://doi.org/10.1038/s41591-018-0057-z

• Kumagai, S., Togashi, Y., Kamada, T. et al. The PD-1 expression balance between effector and regulatory T cells predicts the clinical efficacy of PD-1 blockade therapies. Nat Immunol 21, 1346–1358 (2020). https://doi.org/10.1038/s41590-020-0769-3